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ABSTRACT 

The use of microscopic discrete fluid volumes (i.e., droplets) 
as microreactors for digital microfluidic applications often re- 
quires mixing enhancement and control within droplets. In this 
work, we consider a translating spherical liquid droplet to which 
we impose a time periodic rigid-body rotation which we model 
using the superposition of a Hill vortex and an unsteady rigid 
body rotation. This perturbation in the form of a rotation not only 
creates a three-dimensional chaotic mixing region, which oper- 
ates through the stretching and folding of material lines, but also 
offers the possibility of controlling both the size and the location 
of the mixing. Such a control is achieved by judiciously adjust- 
ing the three parameters that characterize the rotation, i.e., the 
rotation amplitude, frequency and orientation of the rotation. As 
the size of the mixing region is increased, complete mixing within 
the drop is obtained. 
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Affiliated with the CNRS Research Federation FRUMAM (FR 2291). CEA reg- 
istered research laboratory LRC DSM-06-35. 



INTRODUCTION 

Although most microfluidics have been using fluid streams 
as the main means to carry fluids, devices based on individ- 
ual droplets have been proposed as well. In the latter, also 
called "digital microfluidic" systems, "discrete" fluid volumes 
(droplets) rather than continuous streams are used, with the po- 
tential to utilize individual droplets as microreactors | I j. 
Whether fluids are encapsulated within droplets or flow along 
channels, reactions can occur efficiently only in presence of rapid 
mixing. Such mixing conditions are not easy to fulfill due to 
the low Reynolds number of the flow involved, which prevents 
turbulence from taking place. Stirring, in addition to molecular 
diffusion, is needed to stretch and fold fluid elements, thus signif- 
icantly increasing interfacial areas. Whereas some strategies are 
based on complex channel geometry for flows in microchannels, 
active methods (via external forcing) (see, e.g. 1 2 ,3,4,5,6]) have 
also resulted in efficient mixing, especially at very low Reynolds 
numbers 1 7] . The combination of both geometry alteration and 
forcing has been explored as well 1 7, 8,9, 10 , 11 ] . 
Chaotic advection inside a liquid droplet subjected to a forcing 
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has been studied extensively |[T2 13_ 14] 15 16 17 ^ 
tained experimentally by means of oscillatory flows [19 , 20) . In 
this paper as in [21J , we concentrate on controlling both the size 
and the location of the mixing. Indeed, we are interested in the 
ability to control the quantity of fluid mixed, from almost zero 
mixing to complete mixing, as well as in the localization of the 
mixing. On the one hand, researchers have mostly concentrated 
so far on complete mixing as it is often desired for reactions to 
occur uniformly since variations in the mixing could result in 
unacceptable variability in the properties and performance of the 
end-product. On the other hand, incomplete and localized mix- 
ing - when well controlled - could be useful for the synthesis 
of inhomogeneous particles (see, e.g., Janus particles) made of 
two (or more) components whose properties differ. For instance, 
such particles could have a part of their surface hydrophobic and 
some other parts hydrophilic. Such a variation within the same 
particle could result in superior particle properties useful in ap- 
plications such as nanotechnology (self-assembly), biomedicine 
and advanced sensors. 

In this article, we extend the work of (2T\ that focuses on 
unsteady - yet periodic - forcing on a translating droplet and its 
influence on the chaotic regions within the droplet. 
The existence of chaotic behavior in three-dimensional bounded 
steady flows has been shown (e.g. |[T2 13 22^, 23 1). Our work 
is distinct from the latter contributions through the addition of 
unsteadiness, which is crucial to control the chaotic mixing be- 
havior through resonance effects l^i'^J^I • 
The physical system under consideration in this work is de- 
scribed in the first section, and numerical results are exposed in 
the fourth section. We study the chaotic mixing zone by monitor- 
ing its location and its size. Specifically, the creation and control 
of the chaotic mixing regions' location and size are studied qual- 
itatively via Liouvillian sections. A quantification of the size as 
a function of the various parameters involved is also performed, 
thus leading to mixing optimization in parameter space. From 
this quantification of the size, an optimization of the quantity of 
internal mixing is performed. 



PHYSICAL SYSTEM 

Let us assume a Newtonian liquid droplet of sperical shape 
suspended in an incompressible Newtonian fluid whose motion 
consists of a translation and a slow rigid body rotation (see |T3|). 
As in the previous reference, we assume that the interfacial ten- 
sion is sufficiently large for the drop to remain spherical through 
its motion and that the Reynolds number is very small compared 
to 1. Thus, a reasonable approximation is to consider that both 
the internal and external flows are Stokes flows. The boundary 
conditions at the droplet surface can be derived from the conti- 
nuity of velocity and tangential stress conditions. 
The resulting internal flow is a superimposition of a steady 
base flow (non-mixing flow) and an unsteady rigid-body rota- 



tion. Consider a Cartesian coordinate system translating with the 
center-of-mass velocity of the droplet with the orientation of the 
axes such that the unit vector points in the direction of the 
translation and the unit vector Cx lies in the CD — plane. We 
then obtain the following unsteady internal flow: 

= i = zx — £aco(0 sinPj, (1) 
V = j = + 8aco(0 (sin Px - cos P^^) , (2) 
w = z = 1 - 2x^ - 2/ - + 8aco(0 cos Pj, (3) 

In Eqs. ([l][5]), all the lengths and velocities were made dimen- 
sionless by using the droplet radius and the magnitude of the 
translational velocity as the length and velocity scales, respec- 
tively. The rotation is characterized by a maximum amplitude 
8 <C 1, a fixed orientation vector with angle P and a form aco(0 
given by 



<^co(0 = -(1+coscor), 



(4) 



Note that the equations above, together with their derivation, 
are identical to those of 1 13], except that the vorticity vector is 
now time dependent. The time dependency is introduced either 
via the external boundary conditions or through a time dependent 
body force in the momentum equation. 

In practice, this could be realized for instance by creating 
a time-dependent swirl motion in the external flow or by apply- 
ing an electric field that exerts a torque on the drop (see, e.g., 
work on traveling wave dielectrophoresis | 26 | or on electroro- 
tation |j27J). A possible design based on the latter phenomenon 
is proposed in Fig. [T] Electrorotation stands for the spinning of 
an electrically polarized particle (or droplet in our case) while 
the latter is subjected to a "rotating" electric field, that is an ac 
electric field generated by voltages which are out of phase of 
one another. The proposed device would thus consist of a square 
column/channel with electrodes embedded within its four walls 
creating a rotating electric field due to the phase difference be- 
tween the voltages of adjacent electrodes. It is known that such a 
four-pole electrode setting would generate a controlled spinning 
motion, or rigid body rotation, of a drop trapped in the middle 
of the channel. The drop steady translating motion, on another 
hand, could simply be produced by a constant pressure gradient 
along the channel or simply buoyancy in the case of a vertical 
column. In order to vary the orientation of the axis of rotation 
compared to the direction of translation, we propose a series of 
electrodes lying within inclined planes along the length of the 
channel (see Fig. [T]). The angle of the inclined planes with the 
vertical could be adjusted. The periodicity in the angular veloc- 
ity is then obtained by imposing a phase difference to the volt- 
ages applied to the electrodes of adjacent planes, thus generating 
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Figure 2. STREAMLINES IN A CROSS-SECTION OF THE DROPLET 
(WITHOUT ROTATION) AND THEIR FREQUENCIES (\|/) AS GIVEN 
BY Eq. (6). REPRINTED FROM |28] WITH PERMISSION FROM ELSE- 
VIER. 



Figure 1. SKETCH OF A POSSIBLE EXPERIMENTAL APPARATUS 
SHOWING (a) A SQUARE CHANNEL OR COLUMN WITH A SERIES OF 
OUT-OF-PHASE ELECTRODES EMBEDDED IN THE FOUR WALLS. A 
TORQUE ACTS ON THE DROP LOCATED AT THE ORIGIN OF THE 
FRAME OF REFERENCE IN EACH PLANE. SUCH A TORQUE IS IN 
THE DIRECTION PERPENDICULAR TO EACH PLANE, WHICH LIES 
AT AN ANGLE WITH THE STREAMWISE DIRECTION OF THE CHAN- 
NEL (Z). THE VOLTAGES ENERGIZING THE ELECTRODES ALTER- 
NATE BETWEEN THE VARIOUS PLANES (AS SHOWN IN THE FIG- 
URE) SO THAT THE TORQUE SWITCHES BETWEEN POSITIVE AND 
NEGATIVE VALUES AND THE DROPLET UNDERGOES ACCELERAT- 
ING AND DECELERATING ROTATION AS IT TRANSLATES, (b) ACCEL- 
ERATING ROTATION (POSITIVE TORQUE X > 0). (C) DECELERATING 
ROTATION (NEGATIVE TORQUE X < 0). AND CO^ ARE THE VOLT- 
AGE AND FREQUENCY APPLIED TO THE ELECTRODES. THE DASH 
LINE REPRESENTS THE AXIS OF ROTATION. 



acceleration and deceleration phases in the imposed rigid body 
rotation. 



NON-MIXING CASE 

The non-mixing (unperturbed) case or base flow, i.e. 8 = 0, 
is characterized by streamlines that are joint lines of constant 
streamfunction \|/ and azimuthal angle denoted by r^i/^c^. 
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with y(\|/) = vT^^8\j/ and K is the complete elliptic function of 
the first kind. The trajectories of the system ([T][3]) are periodic 
orbits (with period 27r/^l(\i/)) which come in families (labeled as 
Xj/). No chaotic mixing occurs since there is no exponential diver- 
gence in the bulk of the droplet. However, due to the difference 
in frequency a very small mixing occurs in \|/ (but not in ^ due to 
the degeneracy). 



MIXING CASE 

In this work we follow the approach of 1*211 that enables 
the generation of a three-dimensional chaotic mixing region, for 
which we are able to control both the location and the size. The 
method consists in bringing a specific family of the unperturbed 
tori {Tvj;^}^ (composed by a family of periodic orbits along 
the azimuthal angle) into resonance with the periodic perturba- 
tion aco(0 by selecting the frequency CO so that it satisfies the 
resonance condition: 



nO. (\|/„) - CO = 0, for ^ G N. 



(7) 



V = - (x^ +/) (1 -x^ -/ -z^) , ^ = arctan(3;/x) , (5) 



where \|/ G [0, 1/8] and G [0, 27l[. Besides the heteroclinic orbits 
connecting the two hyperbolic fixed points located at the poles 
of the sphere, all other streamlines are closed curves that con- 
verge toward a circle of degenerate elliptic fixed point {x^ -\-y'^ = 
1 /2, z = 0) as \|/ is increased toward the value 1/8. The frequency 
of motion on the streamline r^^(^ is independent of (|) and given 



The control is realized by adjusting the three parameters that 
characterize the periodic rigid body rotation, specifically its the 
maximum amplitude 8, its frequency CO and the orientation of 
the axis of rotation |3 (see Eqs. [l][3]). The amplitude satisfies 
< 8 << 1, where the lower and upper limits correspond to the 
absence of mixing and complete mixing, respectively. 
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Figure 3. LIOUVILLIAN SECTIONS FOR THE FREQUENCIES CO = 
0.95, 1.1, 1.25, 1.40 (a-d), THE AMPLITUDE 8 = 0.05 AND FOR ORI- 
ENTATION a = n/4. THE (RED) LINE INSIDE THE CMRi IS THE 
TORUS r^-i(a))- REPRINTED FROM g WITH PERMISSION FROM 
ELSEVIER. 



Figure 4. LIOUVILLIAN SECTIONS FOR THE FREQUENCIES CO = 
1.25, THE AMPLITUDE 8 = 0.01,0.05,0.10,0.20 (a-d) AND FOR 
THE ORIENTATION a = 7l/4. THE (RED) LINE INSIDE THE CMRi 
IS THE TORUS r^-i(a))- REPRINTED FROM \2Q^ WITH PERMISSION 
FROM ELSEVIER. 



NUMERICAL RESULTS 

The numerical results presented below have been obtained 
by using a standard explicit fourth order Runge-Kutta scheme 

Controlling the location of the chaotic mixing region 

Figures [3j |4] display the Liouvillian sections of the mixing 
system, which consist of two-dimensional projections of time- 
periodic three-dimensional flows by a combination of a strobo- 
scopic map and a plane section (here, the j = plane). In other 
words, the points on the Liouivillian sections are the intersections 
of the trajectories with the plane j = at every period 27l/co. 
The perturbation aco(0 generates two non-negligible three- 
dimensional chaotic mixing regions (see Fig. [3]): one around 
the torus having the frequency CO labeled CMR\ and another one 
around the pole-to-pole axis and near the droplet surface labeled 
CMRn>\. In Fig. |3j we clearly see that the location of CMRi 
varies with the value of CO according to Eq. ([7]). It is important 
to note that CMRn>\ stay anchored around the pole-to-pole axis 
and to the surface due to the flat part of ^^(\|/) close to \|/ = 0. 
This is due to the flat part of 0.{yf) (see Fig. [2]). 
For small values of CO, all resonances are located near the pole-to- 
pole heteroclinic connections (at \|/ = 0, near the z axis and near 
the surface of the droplet, see Fig.jSj. For larger co values, CMR\ 
separates from the chaotic mixing region close to the pole-to- 
pole axis and penetrates deeper into the droplet. In the interval 
< CO < >/2, CMR\ is the largest chaotic region, followed by 
CMRn>\ . As CO is increased further, CMR\ moves toward the 
location of the elliptic fixed point of the unperturbed system by 
following the location of the resonant tori with the frequency CO. 



Controlling and optimizing the size of the chaotic mix- 
ing region 

In this section, we quantify the size of the main chaotic mix- 
ing region, i.e., CMR\ as the three parameters 8, CO and P vary. 
This quantification is performed by computing the the maximum 
variation A\j/ of the stream function \j/ for one trajectory inside 
the CMR\ Whereas the frequency CO of the rigid body rotation 
is mostly responsible for the location of CMR\ (while CMRn>\ 
is always anchored around the heteroclinic orbits), it is the am- 
plitude 8 and the orientation P that mostly determine the size of 
the CMRi. Figure |4] shows qualitatively the size of the chaotic 
mixing regions created by the ^ = 1 and, n> I resonances, i.e., 
CMRi and CMR^^i. In this figure, we see that the size of both the 
CMRi and CMRn>i increases with the amplitude of the pertur- 
bation. It is also interesting to note that around 8 = Emax ^ 0.20, 
the CMRi and CMR^yi join together to cover the entire droplet 
volume. At that point, complete chaotic mixing is obtained. The 
size of CMRi as a function of the frequency CO is depicted in 
Fig. |6] (lower panel). From this figure it is clear that for each 
value of 8 the size reaches a maximum for a certain value co^^^ 
of the forcing frequency. 

The effects of the parameter P on the size of CMRi present 
two distinguishable behaviors one where the size is weakly de- 
pendent and another where the size is strongly dependent. Such 
behaviors are easily observed on the lower panel of Fig. [6] where 
the size of the CMRi as a function of P is displayed. On one hand, 
when P is well inside the two limit values, and 7r/2, the size is 
practically constant. On the other hand, when P gets closer to 
the limits, the size decreases significantly. Such observation are 
qualitatively confirmed in the Liouvillian section in Fig. [5] In 
the first two panels, we see a very small change in size. How- 
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Figure 5. LIOUVILLIAN SECTIONS FOR THE FREQUENCY CO = 
1.25, THE AMPLITUDE 8 = 0.05 AND THE ORIENTATIONS a = 
7r/8,7r/4, 2971/80, 771/16 (a-d). THE (RED) LINE INSIDE THE CM7?i 
IS THE TORUS r^-i(^). 



ever, in the last two panels the size decreases drastically. Having 
such a relation between the size of the CMR\ and the orientation 
of the rotation P can be handy in practice to control the size of 
mixing. In some applications, where the experimenters are faced 
with the dilemma of precisely controlling the size of mixing with 
an imprecise fluctuation of P, setting P well inside the limit value 
and controlling the size through the parameter £ could be the so- 
lution. Then, more restrictive, applications where increasing the 
amplitude of rotation may not be possible (e.g. biomedical ones 
where biological particles need to be handled with care), fixing 8 
and tuning the size of the mixing by varying P around the limit 
value P = 7r/2 could be the solution. One should notice that at the 
critical orientation P = 7r/2 no three-dimensional chaotic mixing 
occurs since \|/ is conserved. 
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Figure 6. UPPER PANEL: NORMALIZED SIZE Axj/ VS. FREQUENCY 
CO FOR AMPLITUDES 8 = 0.01,0.05,0.10,0.20 (a-d) WITH ORIEN- 
TATION P = 71/4; MIDDLE PANEL: NORMALIZED SIZE AXj/ VS. AMPLI- 
TUDE 8 FOR FREQUENCIES CO = 0.55,0.93, 1.28, 1.41 (a-d) WITH 
THE ORIENTATION P = 7l/4; LOWER PANEL: NORMALIZED SIZE AXj/ 
VS. ORIENTATION P FOR AMPLITUDES 8 = 0.01,0.05,0.10,0.20 
(a-d) WITH FREQUENCY CO = 1.25. TOP AND MIDDLE PANELS 
REPRINTED FROM ^ WITH PERMISSION FROM THE AMERICAN 
PHYSICAL SOCIETY 



ACKNOWLEDGMENT 

his article is based upon work partially supported by the NSF 
(grants CTS-0626070 (N.A.), CTS-0626123 (P.S.) and 0400370 
(D.V.)). D.V. is grateful to the RBRF (grant 06-01-00117) and 
to the Donors of the ACS Petroleum Research Fund. C.C. ac- 
knowledges support from Euratom-CEA (contract EUR 344-88- 
1 FUA F) and CNRS (PICS program). 



CONCLUSIONS 

In this work we have shown that chaotic mixing within a 
translating droplet can be obtained by adding a perturbation in 
the form of an oscillatory rigid body rotation. The frequency 
of the latter was selected in order to create resonances with the 
natural frequencies of the system, namely the frequencies of the 
various tori embedded within the droplet. A particularly interest- 
ing feature of the perturbed system lies in the fact that both the 
size and the location of the mixing region can be varied. This 
was achieved by adjusting the control parameters of the pertur- 
bation, specifically the frequency, amplitude and orientation of 
the rigid body rotation. While the frequency is used to target a 
particular location, both the amplitude and the orientation influ- 
ence the size. The latter property introduces additional flexibility 
in the system, particularly for applications which allow for the 
variation of one of these two parameters only. 
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